Pores in Bilayer Membranes of Amphiphilic Molecules: Coarse-Grained Molecular 
Dynamics Simulations Compared with Simple Mesoscopic Models. 
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We investigate pores in fluid membranes by molecular dynamics simulations of an amphiphile-solvent 
mixture, using a molecular coarse-grained model. The amphiphilic membranes self-assemble into a 
lamellar stack of amphiphilic bilayers separated by solvent layers. We focus on the particular case of 
tensionless membranes, in which pores spontaneously appear because of thermal fluctuations. Their 
spatial distribution is similar to that of a random set of repulsive hard discs. The size and shape 
distribution of individual pores can be described satisfactorily by a simple mesoscopic model, which 
accounts only for a pore independent core energy and a line tension penalty at the pore edges. In 
particular, the pores are not circular; their shapes are fractal and have the same characteristics as 
those of two dimensional ring polymers. Finally, we study the size-fluctuation dynamics of the pores, 
and compare the time evolution of their contour length to a random walk in a linear potential. 



I. INTRODUCTION 

Fluid lipid bilayers are the basic material of biolog- 
ical membranes. Pores in such bilayers play an im- 
portant role in the diffusion of small molecules across 
biomembranesiSi^In the last decade, the interest in pore 
formation in bilayer membranes has greatly increased 
with the development of electroporation - in this tech- 
nique, an intense electric field is applied for a short 
time allowing bulky hydrophilic molecules to permeate 
through the lipid membranes of a cell or a vesicle. Ad- 
ditionally, the formation of pores in bilayers is supposed 
to be one key step of the fusion of membranes4*^ 

One difficulty in studying the pores experimentally 
is that they are usually not visible in optical mi- 
croscopy. Yet, they have been investigated using indi- 
rect techniques such as permeation measurements with 
micropipettes ( e.g. Ref. 0) or small angle neutron 
scattering]^ The mechanisms of permeation through a 
membrane obviously depend on the size and the life-time 
distributions of the pores. Such data are difficult to ob- 
tain experimentally because of the small life-time and the 
small dimensions of the pores (less than a millisecond and 
a few nanometers) j^iiS Therefore, numerical calculations 
- either with molecular modelsii^i^iiiiiii^ or with den- 
sity functional theoriesiSiii - have proven useful to study 
the local structure of defects in amphiphilic bilayers and 
lamellar phases. 

Models for pores dynamic3i*i24i2i2& are often based on 
Glaser's model of pore formationi2i He distinguishes be- 
tween three stages. In the first stage, the bilayer thick- 
ness fluctuates because of thermal fluctuations. Some 
hydrophobic tails happen to be exposed to the solvent. 
In the second stage, the solvent spans through the hy- 
drophobic layer, creating a hydrophobic pore. If this pore 
expands, it becomes energetically favorable for the am- 



phiphiles in the edge of the pore to reorient. Finally, in 
the third stage, the pore possesses a hydrophilic edge (see 
Fig. The final hydrophilic pores carry an energy E 



z 
A 



-► r 



FIG. 1; Schematic cartoon of a hydrophilic pore in an am- 
phiphilic bilayer (inspired by Ref. l2ll) . The dark discs are the 
hydrophilic heads of the amphiphiles, the thin lines their hy- 
drophobic tails. The solvent above/under the bilayer and in 
the pore is not represented. In the grey region, the structure 
of the bilayer is perturbed. 



that depends on the contour length c and the area a of 
the pore. Lister— suggested that 



E = Eq + Xc — 7a, 



(1) 



where Eq is equivalent to a core energy, A is the line 
tension of the pore edge, and 7 the surface tension of 
the bilayer. In most cases, A > and 7 > 0. Accord- 
ing to Lister's model [Eq.(l)], for strictly positive surface 
tension, pores larger than a critical size are stable and 
the membrane may break. For tensionless membranes 
(7 = 0, the case we investigate here) the energy increases 
monotonously when the pore grows. Pores of any size are 
unstable, and the membrane is stable. 

The aim of this article is to compare the results of 
molecular dynamics simulations to simple mesoscopic 
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models of pores in bilayers based on Eq. The com- 
parison has turned out to be fruitful for various models, 
and permits to bridge the gap between the descriptions at 
different length and time scales. Miiller and coworkers!^ 
could show that pore formation in stretched bilayers com- 
posed of diblock-copolymers is effectively associated with 
a rearrangement of the amphiphilic molecules in the edge, 
as Glaser suggested. Recently, Talanquer ct aliii con- 
firmed Eq. using a density functional theory for am- 
phiphilic bilayers. 

Here we present a large scale study of spontaneously 
formed pores in a stack of tensionless bilayers. In 
stretched bilayers (7 7^ 0) the pore energy mainly de- 
pends on the pore area, but at zero surface tension 
(7 — 0), the relevant variable is the contour length c 
[see Eq. In particular, we study systematically the 
shapes of pores, their spatial distribution within a bi- 
layer, and their dynamical evolution. 

The paper is organized as follows: in Sec. El we de- 
scribe briefly the model of amphiphilic molecules and the 
molecular dynamics simulation methods which permit us 
to describe the lamellar phas e in th e NPT^ = ensem- 
ble (for more details see Refs. 12312^ . Then, the analysis 
of the molecular dynamics configurations is described. 
The results presented in Sec. liiil are divided into four 
parts. First, in Sec. lIII we present the density profiles 
around the pores. These show that we study hydrophilic 
pores, i.e. pores with a configurational rearrangement in 
the edge. In Sec. IIII Bl the pore positions within a mem- 
brane are investigated. We conclude that the pores do 
not interact unless close to each other. In Sec. IIII CI 
we discuss the size and shape distribution of the pores 
and estimate the line tension A associated to the config- 
urational rearrangement in the pore edge. Finally, we 
trace the time evolution of individual pores. The results 
are described in Sec. IIII Dl and interpreted with a sim- 
ple stochastic modeliifi. We summarize and conclude in 
Sec.HVI 



A. Coarse-Grained Molecular Model 

The model was presented in detail in a previ- 
ous publication^^ and is based on well known similar 
modelSf^SiSi therefore we recall here only its essential fea- 
tures. The bilayers are formed by amphiphilic molecules 
(molecules with a hydrophilic head-group and one hy- 
drophobic tails) in a binary solution with solvent. All 
molecules are represented by one or several soft beads 
(for simplicity, all beads are taken to have the same mass 
to). The solvent is represented by single soft spheres 
(type s). The amphiphilic molecules are linear tetramers 
composed of two solvophobic beads (or "tail beads" , de- 
noted t) and two solvophilic beads (or "head beads" , de- 
noted h). The soft spheres of the amphiphilic are 
connected by bonds. 

Bonded beads are connected by a spring potential2& 
independent of the bead-types: 
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if r < 
if Tfc < r 



This potential comprises a Lennard- Jones type soft re- 
pulsive part dominating for r 0, and an attractive 
part diverging for r —f ri, (the bond length is therefore 
confined between and rf,). The length parameter a is 
our unit of length and the energetic parameter e our unit 
of energy. The bond parameters were fixed at rj, = 2.0 a 
and K = 7.0 e. 

Non-bonded beads interact with short ranged poten- 
tials of the general form 
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(7) -(7) +i 
I [cos(ar2 + /3) - 1] 




if r < 2I/V 

if 2I/V <r<rc 
if Vr < r 



II. MODEL AND METHODS 



At the time-scale available to all-atoms molecular sim- 
ulations (about 100 ns). the spontaneous formation of 
a pore in a lipid bilayer can be considered as a rare 
event. In fact, the formation of a pore in a DPPC bi- 
layer was studied a few years ago with all-atom molecu- 
lar dynamicsi^ of 256 lipids during 0.5 fis, but this pore 
was observed in a far-from-equilibrium situation with a 
relatively small system. To explore the properties of 
thermally activated pores in equilibrium bilayers, we use 
a coarse-grained molecular model which describes rela- 
tively thin bilayers. In such simulations, many pores ap- 
pear spontaneously during one simulation run. 



This potential comprises again a Lennard- Jones type soft 
repulsive part, and a short-ranged attractive part. The 
range of the potential is fixed at 1.5 cr. The pa- 
rameters a and f3 are fixed such that potentials and 
forces are continuous everywhere {a = 7r/r^ — 2^/'^a^ and 
(3 ~ 27T ~ r'^a) . The energetic parameter (p determines 
the depth of the potential; it depends on the types of 
the interacting beads. For pairs that "dislike" each other 
(ts and th), we choose cj) — and so that the interac- 
tion is purely repulsive. For the sake of simplicity, for all 
the other pairs {ss, sh, hh, and tt), the potential depth 
(j> of the pair interactions is the same {(p = 1.1 e). For 
a fixed e, the self-assembly is driven by the increase of 
(j) only. The choice of (j) = 1 . 1 e ensures that we simu- 
late the liquid crystalline lamellar phase of the tetrameric 
amphiphilesi^ 
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In the following, lengths shall be given in units of a, 
energies in units of e and masses in units of m. This 
gives the time unit r = (mcr^/e)^/^. Typical orders of 
magnitude of our units are for example e ~ 5 • 10^^^ J, 
m 10~^^ kg, cr ~ SA, and r 10"^^ s. In this 
case, one solvent bead represents roughly three water 
molecules, and the tetrameric amphiphiles typically rep- 
resent a symmetric neutral surfactant like CqEO^. Such 
surfactants are much smaller than usual biological lipids. 
As a consequence, bilayers diluted in a large amount of 
solvent are not stable in our model (they separate into 
micelles). Thus, a quantitative comparison of our sys- 
tem with realistic lipid biomembranes is not our goal: 
we focus on the general properties of pores in bilayer 
membranes. 

To study stable bilayer membranes, we simulated self- 
assembled La lamellar structures: a stack of amphiphile 
bilayers parallel to each other, separated by layers of 
solvent (sec Fig. At the temperature of the simu- 
lation, these bilayers arc two dimensional fluidsi2i24 In 
the present work, a smectic phase composed of fifteen 
bilayers containing several thousands of molecules each 
was simulated over 10^ t. 




FIG. 2: Snapshot of the system (30 720 /i2t2 and 30 720 
solvent beads simulated in the NPT ensemble with P/e — 
2.9 cr"^, ksT = e, and c^/e = 1.1). The dark beads are 
solvophobic (type t), the light beads are solvophilic beads or 
solvent beads (type h or s). 



B. Simulation Details 

We have analyzed the same simulation runs as in a 
previous related article;^^ in which the lamellar La phase 
and its elastic properties are characterized in more detail. 
Here, we focus the analysis on the defects appearing in 
the bilayers due to thermal fluctuations. We have stud- 
ied the model in the (iVP„Ptr)-ensemble (constant num- 
ber of particles, constant pressure normal and tangential 
to the bilayers, and constant temperature) with molecu- 
lar dynamics simulations. (N): The lamellar phase was 
studied at an amphiphile fraction of 80% of the beads 
(one solvent bead per /i2i2)- The system contained 30 720 
tetramers and 30 720 solvent beads, which formed fifteen 
parallel bilayers of about two thousand molecules each 
(see Fig.EJ. 

(P): The normal and tangential pressure components 
Pn and Pt were kept constant using the extended Hamil- 
tonian method of Anderseni^Sj^S The box shape is con- 
strained to remain a rectangular parallelepiped. The box 
dimension perpendicular to the bilayer (L^) and tan- 
gential to the bilayers {Lx,Ly) are coupled to two sep- 
arated pistons. The thermal- averaged box dimensions 
were L^ = Ly = 43.4 ± 0.1 cr and L,_ = 31.9 ± 0.1 cr. We 
imposed the two pressure components separately rather 
than the total pressure P = (P„ + 2Pt)/3 because of 
technical reasons: the mechanical equilibrium is reached 
earlier, the orientation of the bilayers is stabilized, and 
the surface tension is controlled. Since we studied a 
bulk lamellar phase, we imposed an isotropic pressure 
{Pn ^ Pt ~ P)- More details on the simulation algorithm 
and its parameters are given in Ref.l23i (T): The temper- 
ature was controlled by a stochastic Langevin thermostat 
that has been described earlier and applied to very sim- 
ilar modclSiSSiSLffi 

The dimensionless pressure and temperature were fixed 
at Pcr^/e = 2.9 and ksT/e = 1.0. For the chosen dimen- 
sionless potential depth, 0/e = 1.1, the density fluctuates 
around 0.85 beads per unit volume. With these typical 
parameters, lamellae form and order spontaneously. But 
this process requires at least 30 000 r. We have there- 
fore imposed the orientation of the lamellae in the initial 
configurations. They were constructed so, that fifteen 
bilayers separated by solvent layers were stacked in the 
z-direction. These configurations were then relaxed for 
10 000 T. During that time, the interlamellar distance 
adjusted to its equilibrium value, the shape of the flexi- 
ble box changed accordingly, but the director remained 
basically aligned with the z-direction. 

Data were then collected over 10^ r. We veriflcd that 
the pressure tensor obtained at equilibrium was diag- 
onal and isotropic, and that the surface tension 7 ^ 
{Lz{Pn — Pt)) was negligiblejSi For the time indepen- 
dent analysis (Sees. IIII Al ITiTbI and ITlI C|) . we have used 
400 configurations separated by 250 r. For the time de- 
pendent analysis (Sec. IIII Dp . we have used four series 
of 400 configurations separated by 1 t. These series were 
separated by 10 000 r. 
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C. Pore Detection and Analysis 

Generally, several types of point defects exist in smcc- 
tic phases. A highly aligned La phases may contain pores 
in the bilayers, but also necks or passages between two 
bilayersi^i*^ Any of those defects changes the topology 
of the lamellar phase: pores in the bilayers connect two 
neighboring solvent layers, fusion between neighboring 
bilayers connect the bilayer matrices, and passages con- 
nect both neighboring solvent layers and neighboring bi- 
layers. One possibihty to detect a defect, and to distin- 
guish between these defects is to analyze the topology 
of the lamellar phasei^i^^ Using a cluster-algorithm on 
the simulated lamellar phase^^ we observed that in our 
model La phase, necks and passages appear very rarely 
(in less than 1% of the configurations). So we can focus 
on the pores in the bilayers. 

Figure O shows a representative snapshot of the tail- 
beads of one membrane. On the snapshot, five pores are 




FIG. 3: Snapshot of the head-beads of one bilayer (top view). 



clearly present. 

For a systematic analysis, we defined a pore as a region 
of the inner part of the bilayer were relatively few tail- 
beads are present (i.e. the relative density of tail beads 
is less than a given threshold) . The nth bilayer is defined 
by its position hn{x,y) and its thickness tn{x,y) (Monge 
representation)^ In practice, only discrete values of x 
and y were considered {x = n,j.Lr^/N,j. and y = UyLy/Ny 
where nx,y are integer values, — Ny = 32, and 
Lx = Ly ^ 43 a). Therefore, all the obscrvablcs are mea- 
sured on a two dimensional grid of 32 x 32 "plaqucttes" 
of area dx x dy where dx = Lx/Nx and dy = Ly/Ny. As 
the dimensions of the simulation box fluctuate, the mesh 
size also fluctuates: dx = dy ^ 1.3 ± 0.05 cr. In the fol- 
lowing, the notation "plaquette" indicates that area are 
expressed in units of the plaquette area and the lengths 
in units of the mesh size. 

The analysis is divided into two steps (see Appendix 
^for more details): First, we have determined the local 



positions and thickness of the membranes in every config- 
uration from the relative densities of solvophobic beads 
Ptaii{x, y, z). For each point {x, y), the relative density of 
tail beads ptaii{x,y, z) oscillates as a function of z, with 
one maximum per bilayers. The position hn{x,y) and 
thickness i„(a;,y) of a given membrane were determined 
as the mean and the difference of the two z values where 
Ptaiiix,y, z) equals a threshold value (po ~ 0.7) around 
the appropriate maximum. 

Second, the positions {x,y) where the thickness of the 
membrane n is zero or undefined are attributed to pores. 
We call them "pore-positions" . For each membrane, the 
pore-positions {xi,yi}PP are assembled into two dimen- 
sional clusters, the "pore-clusters". Two pore-positions 
belong to the same pore-cluster if they share at least one 
vortex of the two-dimensional Nx x Ny grid. 




FIG. 4: Result of the pore analysis on the membrane depicted 
in Fig. |S| The single isolated square in the bottom left of the 
image is of the area of a single plaquette. 

The result of such an analysis is represented on Fig. 01 
In this figure, the dots represent the centers of the tail- 
beads belonging to the membrane (seen in Fig.OJ. Each 
pore-position is represented by a grey square. The re- 
sulting grey patches are the pore-clusters. Each pore 
visible in Fig. |2| appears under the form of a pore-cluster. 
But additional small pores-clusters that were not visi- 
ble on the snapshot also appear in the analysis. It seems 
that some tail-beads (black dots) are still present in those 
pores. In practice, it was difficult to decide whether these 
small pore-clusters are noise, fluctuations of the bilayer 
thickness, or hydrophobic pores. In particular, no mini- 
mum life-time was found for these pores of minimum size. 
Therefore, the pore-clusters composed of one plaquette 
only were disregarded in the following analysis (unless 
noted otherwise). 

For each of the pore-clusters, which we call only 
"pores" for simplicity, the mean position rem, the ma- 
trix of gyration g"^ , the area a and the contour length c 
were calculated. Again, these observables are defined for 
clusters of pixels on the two dimensional Nx x Ny grid. 
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The mean position is calculated via 



TAILS 



(4) 



where the sum runs over the Upp pore-positions of the 
pore cluster, and is the position of the pore-position 
number i. 

The gyration matrix is given by 

-1 ^pp 

g"^ = .90 + V(r,; - r,„,)"(r, - r™,)^ (5) 

7? ^ ^ 



where a and P denotes axes of the grid, and go the gyra- 
tion matrix of one plaquette, which we approximated by 
one fourth the identity matrix (in plaquette area). 



III. SIMULATION RESULTS 

In each bilayer of area L^Ly ^ 1850 cr^, the average 
number of pores is 9.9 ± 1.3. Their total area corre- 
sponds to approximatively 1% of the projected area of 
the bilayer. Among these ten pores, 5.7 ± 1.0 pores on 
average have the minimum size of one plaquette (with an 
area of about 1.7 a^). 



A. Composition Profiles through the Pores 

Figure |31 shows the composition in head, tail and sol- 
vent beads around the center of relatively small pores (of 
area two to four plaquettes) . They are represented in the 
cylindrical coordinates with the center of the pore as ori- 
gin: ptatiir, z), pheadir, z) , Psoiventi^j z) , whcrc thc z-axis 
is perpendicular to the bilayers (see Fig.P). We averaged 
over the directions of r in the plane of the bilayer and over 
the two sides of the bilayers for the values +z and —z, 
so the composition is represented as a function of |r| and 
\z\. 

Far from the pore centers, at distances r > 4tT, the 
typical bilayer structure is found. The tails are segre- 
gated inside the bilayer (z < 2 a). For intermediate z 
(around 2a), the head-beads shield the tail-beads from 
the solvent. The solvent is concentrated near the mid- 
bilayer plane {z = 3.2 cr). As defined by our detection 
algorithm, the pore (r ^ 0, z < 3ct) is characterized by 
the absence of tail-beads inside the bilayer. The head and 
solvent distributions show that the interior of the pore is 
occupied by head-beads rather than solvent-beads. Even 
inside the pore, the head-beads shield the tail-beads from 
the solvent. We conclude that pores with an area a larger 
than two plaquettes are of "hydrophilic type". As ex- 
pected, the same observation is made for larger pores 
(data not shown). To summarize, the amphiphiles of the 
pore edge are reoriented in our model, even for the small 
pores (a = 2,3,4 plaquettes). Notably, such a reorienta- 
tion was also observed by Miiller and co-workers in the 
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FIG. 5: Composition profiles around the centers of pores 
with areas between two and four plaquettes (length unit: a). 
Note that the inter-bilayer distance is 6.4 cr, so the mid-bilayer 
plane (between a bilayer and its next neighbor) is at z = 3.2 o 



pore edges of amphiphilic bilayer of block-copolymers4^ 
As a consequence, it is reasonable to use Eq. ^ to inter- 
pret our simulation data: the energy of a pore is expected 
to increase with the contour length of the pore. 



B. Positions Correlations of the Pores 

About four pores of area larger or equal to two pla- 
quettes are present in each simulated membrane of area 
Lj. X Lj, ~ 43 X 43(7^. Do these pores interact? A priori, 
they are at least subject to hard core repulsion, because 
two pores cannot overlap (overlapping pores would be 
considered as a single pore). In this section, we try to 
detect whether there exist additional soft interactions be- 
tween the pores. 
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FIG. 7: Principle of the Minkowski analysis (see text for more 
explanations) . 
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FIG. 6: Spatial pair correlation function of the center of the 
pores. As we neglected the correlations between the distribu- 
tions of pores in different membranes of the same configura- 
tion, the errorbars are under-estimated. 



The most straightforward approach to this problem is 
to look for spatial pair correlations between the pores 
(see Fig-EJ. In this analysis, we have taken into account 
all the pores, even the smallest ones (a = 1 plaquette). 
Despite the noise of the data, two tendencies are clear. 
At large distances, say about r > 10 fi, no correlation 
is observed. At small distances (r < 5 a), the short- 
range repulsion between the pores is reflected by a de- 
pletion in the pair correlation function. At intermediate 
distances (r ~ 8 c), the data suggest a slight positive cor- 
relation, i.e. the probability that pores open up seems to 
be slightly increased in the vicinity of other pores. The 
statistics is however not sufficient to justify a conclusive 
statement. 

An alternative, complementary method of analyzing 
spatial distributions of patterns in physical systems has 
been proposed by Mecke and other Q?'^i?^i?^i?Ti??i?? : the 
evaluation of Minkowski functionals. Minkowski func- 
tionals have been used in various areas of mathematics, 
chemistry and physics to analyze high-order correlations 
of spatial distributions. The present analysis is similar 
to the one developed and applied by Mecke, Jacobs and 
co-workers^Mli to study the spatial distribution of de- 
fects in a thin film of polymers. The analysis can be 
decomposed into two steps (see Fig. First, on each 
of the N points of the distribution, a square of size / is 



fixed (the gcrm)^ These N germs form a two dimen- 
sional image, the pattern. In two dimensions, the three 
dimensionless Minkowski functionals too,i.2 of a pattern 
are proportional to its surface A, contour length C and 
Euler characteristic tuq = A/ At where At the to- 
tal area of the surface, nii ~ C/\/A[N, where N is the 
number of points in the distribution, and m2 = x/-^ • 
We used an algorithm described by MilchielsenSi to cal- 
culate the Minkowski functionals of a digitalized pattern 
(without periodic boundary conditions). Finally, the di- 
mensionless Minkowski functionals arc represented as a 
function of the dimensionless coverage factor x = PN/At 
where P is the area of the square germ. 

Figure |H1 shows the three Minkowski functionals ob- 
tained from about 1000 pore distributions containing ex- 
actly six pores of area larger than two plaquettes (sym- 
bols). We compare these results with three reference dis- 
tributions (curves). 

The first reference is a random distribution of fixed 
germs on an infinite surface (Poisson), for which the 
Minkowski functionals arc known analytically^: 

mQ^x) = 1 — cxp(— x) (6a) 
TOi(x) ~ 4v^exp(— x) (6b) 

TTl2ix) — (1 — x) CXp(— 2") (6c) 

The pore distributions obviously differ from a Poisson 
distribution on an infinite surface. We attribute the dis- 
crepancy to two effects: finite size effects, and repul- 
sion effects. The former are due to the finite size of our 
quadratic grid, and to the finite area of the membranes. 
Explicit expressions for the Minkowsky functionals in a fi- 
nite "observation window" on an infinite surface are given 
in Ref. 13 Here, we have additional boundary effects be- 
cause only germs inside the window contribute to the 
pattern. Therefore, we have calculated the expected fi- 
nite size effects numerically. 

Two reference distributions on a finite grid of 32 x 32 
pixels were considered. In the first one, which we call "fi- 
nite random" (FR) distribution, a fixed number of points 
was distributed randomly. The second one, denoted "fi- 
nite self-avoiding" (FSA) takes into account the hard- 
core repulsion between pores, i.e. pores may not touch 
or overlap (see Appendix). As Figure |S1 shows, the differ- 
ences between FR and FSA distributions are relatively 
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small, but observable. We present here the results for bi- 
layers containing six pores, but the same conclusions are 
drawn for different densities (the data are not shown; as 
expected, the influence of self-avoidance increases with 
the density of points). 
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FIG. 8: Dimensionless Minkowski functional as a function of 
the dimensionless coverage factor (6 points on a 32 x 32 grid) . 
"Poisson" corresponds to Eq. |S| FR to the numerical results 
obtained with 1000 Finite Random distributions, and FSA to 
the numerical results obtained with 1000 Finite Self- Avoiding 
distributions (see text and Appendix). 



The results obtained with the pore distributions of the 
simulations are closer to those obtained with the FSA 
distributions. Therefore, both finite size effects and self- 
avoidance are observable in the Minkowski analysis, and 
these two effects seems to be sufhcient to account for the 
data. 

To conclude, the pores are not distributed randomly in 
the membrane. The spatial distribution of pores is com- 
patible with a simple model, in which the pores interact 
only through a hard-core repulsion. Additional soft in- 
teractions between pores may be present, but these do 
not influence the overall pore distribution significantly. 

We have also studied the correlations between pores in 
different membranes. Not surprisingly, a small correla- 
tion is observed, i.e. the probability that a pore opens 
up on top of another pore is slightly increased. But we 
did not explore this effect throughout. 



C. Size and Shape of the Pores 

In the following the area of the pores, their contour 
length and their radius of gyration are analyzed. In this 
section, all pores have been taken into account, even the 
smallest ones [a — \ plaquette). 

As shown in Section [111 Bl the spatial distribution of 
the pore within a membrane were in good agreement with 
the hypothesis that the pore interact only via a hard- 
core repulsion. In the following, we consider the pores as 
independent. Within this hypothesis, the contour length 
distribution P(c) of the pores can be compared to the 
Bolzmann factor 5(c) exp[— £'(c)/(fcBT)], where E{c) is 
the energy to create a single pore of contour length c 
and g{c) is the degeneracy. We calculated g(c) for the 
particular quadratic grid of our analysis (see Table QJ. 

The ratio P{c)/g{c) is shown in Fig. |51 in a linear- 
log plot. The approximated energy of pore formation 
— In [P(c)/g(c)] is well described by a linear function. 
Fitting the model E{c) = E'o -I- Ac to the data yields the 
line tension A = 0.7±0.1 ksT-a^^ . This value is approx- 
imately 5 • 10~^^ J • m~^, which agrees reasonably with 
the values calculated by May^ for the excess free energy 
of the packing rearrangement of amphiphiles in the edge. 
Previous resultaSiiSii^ report line tensions in the range of 
10^^^ J • m^^ to 3 • 10~^^ J • m^^, i.e. larger than the 
present value. In these references, the line tension in- 
cludes also the excess free energy necessary to transfer 
the amphiphiles from a reservoir to the edge of pore (the 
chemical potential of the amphiphiles in the grey region 
of Fig. This contribution (typically 10""'^^ J • to~^) is 
proportional to the surface tension of the bilayeriiSiiS and 
vanishes in the present case. 

In brief, the size distributions computed with the sim- 
ulations are compatible with the usual mesoscopic model 
of pores energetics and permit to compute the approxi- 
mate line tension of the pore edge. 

The shape of the pores have been studied with the 
two-dimensional gyration matrix of the pore-clusters [see 
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14 


16 


g(c) 
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9 36 


168 


715 


2000 



TABLE I: Degeneracy of the contour length c of clusters of 
pixels connected by at least one vortex. As explained in Sec- 
tion 8- pore is defined as a cluster of pixels connected by 
at least one vortex. In our calculations, the directions —x, x, 
—y and y are distinguishable. We neglect the finite size of the 
grid (a square of 32 plaquettes), therefore all positions of the 
pores on the grid are considered eis equivalent. Because of the 
high cost of such calculations, the value are exact only up to 
c = 14. The value for c = 16 is an under-estimate. 
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FIG. 9: Probability distribution function of the contour of the 
pores P(c), divided by the degeneracy of the contour length 
g{c) (linear- log plot). 
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FIG. 10: Square of the radius of gyration of the pores, as a 
function of their area. 
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Eq. Q]. The two (positive) eigenvalues of the gyra- 
tion matrix are denoted pi and p^. The sum of the 
eigenvalues Rg = Pi + pi is the square of the ra- 
dius of gyration of the pore, and the relative difference 
I Pi ~ PilZ-^g) its asymmetry - it is zero for a circular 
pore and tends towards one when the pore is elongated 
in one direction. Figure [TUI represents the radius of gy- 
ration of the pores as a function of their area. As ex- 
pected, the square of the radius of gyration increases lin- 
early with the area of the pores. A least-square linear 
fit yields - 0.17a + 0.29 (both in units of cr^). The 
proportionality factor (0.17) is slightly larger than the 
proportionality factor obtained for a homogeneous disc 
([27r]-i ~ 0.15). 

As illustrated in Fig. ^] the asymmetry of the pores 
does not vary significantly with the size of the pores, ex- 
cept for the smallest pores, whose asymmetry is imposed 
by the finite mesh size of the grid. 

This suggests that only one length-scale is sufficient to 
describe the pore dimension. Notably, the average asym- 
metry is not zero (0.35 ±0.05). In fact, since the bilaycrs 
are flaccid, there is no reason why the pores should be 
circular. Miiller and co-worker9& emphasized the impor- 
tance of the asymmetric shape of pores on the mechanism 
of fusion between two parallel bilayer membranes. 

The pore shape can also be studied via the corrcla- 
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FIG. 11: Asymmetry of the pores as function of their area. 

tion between the area and the contour length. Figure [T^ 
shows the pore area a as a function of the contour length 
c in a log-log representation. The simulation data are 
well represented by the equation 

a oc c3/2. (7) 

Such a relationship is consistent with the results obtained 
by Shillcock et ali^ for their study of the proliferation 
of pores in fluid membranes. As Shillcock et al. remark, 
two dimensional flaccid vesicles^ and two-dimensional 
ring polymers^ show the same behavior as pores in flac- 
cid membranes. These different objects may be modeled 
as closed, self- avoiding; planar random-walks whose en- 
ergy depends only on the number of steps. As a conse- 
quence, the correlation between the area a and contour 
length c observed in our simulations [see Eq. Q] is con- 
sistent with the simple model of the pore energy given 
by Eq. ^ (with 7 = 0). We emphasize that the present 
study and Ref. are not using the same types of model. 
Shillcock et al. simulated pores in membranes at a mcso- 
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FIG. 12: Area of the pores as a function of their contour 
length (log- log plot). 



scopic length-scale, with Monte Carlo simulations based 
on Eq. whereas in our simulations, Eq. |^ emerges 
from the molecular model. 



D. Dynamics of size fluctuations 

In the previous section, we have shown that static 
properties of the pores (size and shape distributions) are 
in very good accordance with Eq. (1). In the following 
section, we investigate the dynamics of the pore size by 
comparing the simulation results to a simple stochastic 
model based on Eq.(l). First we determine the time-scale 
of the pore dynamics. Second, we present the stochastic 
model and we check the relevance of the assumptions 
of the model for our simulation results. Finally, the 
"shrinking-time" distribution of the pores of the model 
is fitted to the simulation results and the resulting pa- 
rameters are briefly discussed. 

To characterize the short-time dynamics of the contour 
length c, we have measured the number of jumps ri(c, dc) 
from c to c + dc during one time interval (At). The 
probability P{c, dc), that a jump with the initial contour 
length c has the amplitude dc is then defined as 



P(c, dc) 



ri(c, dc) 

Edc^(c,C^c)' 



(8) 



An appropriate numerical value for At is slightly lower 
than the typical time that a pore needs to change its size. 
The time-scale of pore dynamics is expected to be of the 
same order of magnitude as that of configurational rear- 
rangements of the amphiphiles, i. e., 10^^^ s to 10^^^ s 
(close to It). To estimate this time-scale, we studied the 
correlation time of average quantities like the total con- 
tour length, and the total area of pore per bilayer (see 
Fig. ll3|l . As expected, the correlation times range around 
2 T. Given this first estimate of the time-scale dynamics. 




FIG. 13: Autocorrelation functions of the total areas and 
contour lengths as a function of time (time unit: r). The 
pores of area a = 1 plaquette are not included in the analysis. 



we studied the time-evolution of individual pores every 
It. 

We compare our results with a simple stochastic model 
that has been solved analytically: a random walk in a lin- 
ear potential (RW-LP model) The RW-LP model de- 
scribes a random walk in a semi-infinite one-dimensional 
space with discrete states labeled n = 0,l,2,...,oo. For 
the interpretation of our simulation results, the vari- 
able n{t) represents one half of the contour length c{t). 
The time evolution consist of discrete jumps of ampli- 
tude \dn\ = 1, at an average rate W. The probabil- 
ity to hop towards the upper neighbor {n — > n + 1) is 
p+ = (1 — b)/2. In the other direction (n — > n — 1), it 
is = (1 + 6)/2. The bias < < 1 measures the 
tendency to walk towards n = 0, where the walker dies 
(absorbing boundary) . In our interpretation, the param- 
eter b is proportional the line tension of the pores. For 
this model, the probability Q{n, t\no) that a walker start- 
ing from the state no reaches the state n < no for the first 
time after having walked during the time t iaS^ 



Q(n,t|no) = 



t \ l-b 



-Wt i 



{Wt^/T^) , (9) 



where Ivt{v > 1) is the modified Bessel Function of the 
first kind, and W and h are the two parameters of the 
model. In Eq. the time t is the time needed by the 
variable n to shrink from ng > n to n, so we call it 
"shrinking time" and Q "shrinking-time distribution" . 

Of course, the simulation data are more complicated 
than the RW-LP model, (i) The model supposes that 
within the time W~^, only jumps of amplitude \dc\ = 
2 occur. In the simulations, its is obviously not the 
case . Nevertheless, the jump probability P(c, dc) de- 
creases when the amplitude of the jump \dc\ increases 
(see Fig. HHl. For our choice of At = It, 80% of the 
jumps with dc ^ Q correspond to the amplitude \dc\ = 2. 
One might expect that if the time of observation is small 
enough, we should observe exclusively jumps of the small- 
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FIG. 14: Probability of jumps from c to c + dc as a function of 
dc for possible values of the initial contour length c (linear-log 
scale) . 



est amplitude. We find however that even for the obser- 
vation time At — 0.1 r. some jumps with large amplitude 
\dc\ > 6 appear. The frequency of jumps is thus broadly 
distributed. 

(ii) The model supposes that the effective potential 
(— In P(c)) is a linear function of c. In the simulation, 
it is not linear because of a large cntropic contribution 
to the free energy (data not shown here, see Ref 12^) . 
Nevertheless, the effective potential increases with the 
pore size; In Fig. the energetic bias is clearly seen: 
for a given amplitude \dc\, the probabilities to jump with 
dc < are much larger than for dc > 0. 

(iii) The model does not describe the complex behav- 
ior of small pores, especially the dynamics of hydrophilic 
pore formation. 

After discussing the assumptions of the RW-LP model, 
we show that this stochastic model, yet simple, fit the 
simulation data. We have measured Q{n,t\no) for the 
simulations using no = 5 and n = 4. This choice per- 
mits to restrict ourselves to large pores (c > 8); there- 
fore, we do not describe the dynamics of the formation 
of the pores, but only of their shape and size fluctu- 
ations. As expected from Eq. ©, other choices of 
{n, np} with uq — n = 1 and ng > 5 give very similar 
results. Fig. 1151 shows the comparison between the sim- 
ulation data (symbols) and Eq. @ with the parameters 
b = 0.2 and W = 0.5t~^ (curve). Despite the simplicity 
of the model, the agreement is very good and the fitted 
parameters are reasonable: The value of the frequency 
W = 0.5 is consistent with the previously estimated 
correlation time of the size fluctuation (2 t). The value of 
the bias b = 0.2 is also consistent with the slope s of the 
free energy — In P(c) as a function of c (s = 0.1 fcsT per 
plaquette^ To conclude, not only the statics, but also 
the dynamics of the pore size is in good agreement with 
Lister's model. 
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FIG. 15: Probability that a pore of contour c = 10 reaches the 
contour c = 8 for the first time after having fluctuated during 
the time t (observation time At = Ir), in a linear-log plot. 
The line represents Eq. Q with no = !3,n = 4,W = 0.5 r~^, 
and b — 0.2. We had to renormalize the pore distribution to 
fit it to Eq. because it was normalized on the range t £ 
[2, 450] whereas Eq. Q is normalized on the range t G [0, oo]. 



IV. CONCLUSIONS 

We have investigated a bulk lamellar phase in an am- 
phiphilic system by molecular dynamics simulations, us- 
ing a phenomenological off-lattice model of a binary 
amphiphile-solvent mixture. The system was studied 
in the (NP„PtT)-ensemble using an extended Hamil- 
tonian which ensured that the pressure in the system 
was isotropic. Therefore, the membranes had no surface 
tension. At high amphiphile concentration, (80% bead 
percent of amphiphiles), the amphiphilic molecules self- 
assemble into a lamellar phase, i.e., a stack of bilayers. 
The pores appearing spontaneously in the bilayers were 
studied. 

The molecular structure of the pores shows that the 
amphiphiles situated in the rim of the pore reorient: the 
hydrophilic heads shield the hydrophobic tails from the 
solvent. 

The mean effective free energy of a single pore was 
computed from the distribution of the contour-lengths c 
of the pores. Taking into account the cntropic effect of 
shape fluctuation, we could fit the simple model E(c) — 
Eq + Ac to the simulation data, and estimate the line 
tension A of the pores. 

Without surface tension, the pores are not circular. 
The relationship between the area a of the pores and their 
contour-length c is well described by the law a oc c^^^. 
This relationship was found for other two-dimensional 
objects, whose energy depends only on their contour 
length (models of flaccid vesicles^ and of self-avoiding 
ring-polymers) 4i The analogy between our simulation 
results and these "random-walk" models suggest that in 
our analysis, the "bending energy" of the pore edge is 
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negligible. 

Finally, the shrinking-time distribution Q(ri,t\no) de- 
creases relatively slowly with the time t (it can be fitted, 
for example, by the power law T~^). The long tail of this 
probability distribution indicates the existence of par- 
ticularly stable pores. Presumably, these correspond to 
pores that have grown very large before shrinking again. 
Despite its simplicity, the model of one-dimensional ran- 
dom walk in a linear potential (RW-LP)2^ reproduces 
nicely the distribution of the shrinking times observed 
for the pores. 

To summarize, we have studied in great detail the pore 
statics and dynamics in a molecular model, and found 
that our data are in excellent agreement with the predic- 
tions of a mesoscopic line tension theory E(c) = Eq + Ac. 
Our results establish that such simple phenomenological 
models do indeed describe many aspects of real pores in 
self-assembled membranes. We have demonstrate that 
for a relatively simple, coarse-grained model, which does 
however treat the particles on a microscopic level. Thus 
we believe that our conclusion will also hold for more re- 
alistic models, as long as long-range interactions can be 
neglected. 
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APPENDIX A: ALGORITHM TO DETECT AND 
ANALYZE THE PORES. 

This appendix describes how we determined the local 
positions of membranes in the lamellar stack and the po- 
sition of the pores in the membranes. 

1. The space is divided into N^NyN^ cells of size 
{dx,dy,dz) with = Ny ^ 32 and ~ 92. 
For a density of 0.85 particle per volume unit, 
dx = ~ 1.3 (7 and ~ 1.0 a. The size of cells 
may slightly vary from one configuration to another 
because the dimensions of the box dimensions vary. 

2. The relative density of tail beads in each 
cell is calculated as the ratio ptaii{x,y,z) = 
Ntaii{x,y,z)/Ntot{x,y,z) where Ntaii{x,y, z) is the 
number of tail beads in the cell {x,y,z) and 
Ntot{x,y, z) the total number of particles in this 
same cell. 

3. The membranes are defined as the space where the 
relative density of tail beads is higher than a thresh- 
old {ptaii{x , y , z) > po). The choice of the thresh- 
old depends on the mesh size in x and y directions 
{dx = Lx/Nx and dy = Ly/Ny). Typically, we used 



Po from 0.65 to 0.75 (80 % of the maximum relative 
density of tail beads; values for dx = dy ^ l.Scr). 

4. The cells that belong to membranes are associated 
into three dimensional clusters: Two membrane- 
cells that share at least one vortex are attributed 
to the same membrane-cluster. Each membrane- 
cluster defines a membrane. This algorithm identi- 
fies membranes even if they have holes. In the pres- 
ence of necks between adjacent membranes (local 
fusion), additional steps have to be taken in order 
to find the necks until one membrane-cluster per 
membrane is found (not detailed here). 

5. For each membrane n and each position {x,y), the 
two heights /i™™(a;,y) and h'^°-^{x,y) where the 
density ptaii{x,y, z) equals the threshold po a-re es- 
timated by a linear extrapolation. The mean posi- 
tion and the thickness are then defined by 

hn{x,y) = \ [hr''{x,y)+K"'{x,y)\ 

tn{x, y) = \ [/ir"(x, y) - hT-{x, y)] (Al) 

6. The positions (x, y) where the thickness of the 
membrane n is zero or undefined are considered as 
"pore-positions" . For each membrane, an ensemble 
of pore-positions {xi^yi}^^ is obtained. 

APPENDIX B: FSA DISTRIBUTIONS 

We sketch here the principle of our construction of "fi- 
nite self-avoiding" (FSA) distribution of NPOINTS points 
of size size on a grid of booleans template of dimen- 
sions NX*NY. Each boolean template [i , j] can take the 
value FREE or OCCUPIED. A schematic algorithm can be 
written in the following way (not in a real programming 
language, and without any control about the feasibility 
of the task !): 

function ConstructFSA (NX , NY , NPOINTS , size) 
n = 0; 

Initialize (template , NX , NY , FREE) 
while (n<NPOINTS) 
xn = NX*random() 
yn = NY*random() 
if (template [xn,yn] == FREE) 
x [n] = xn 
y[n] = yn 

Set_occupied (template , xn , yn , size) 
n = n+1 

end_if 
end_wliile 
end_f unction 

where random () is a random number generator whose 
output ranges between and 1. At the end of the loop 
(it is ever reached !), the vectors x and y contains the 
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coordinates of the points of the distribution. The func- 
tion Set_occupied(template,i, j ,size) attributes the 
label OCCUPIED to all the sites of the template around 
the site (i , j ) . 

function Set_occupied(template,i, j ,size) 
a = integer ( (size+l)/2) ; 
for(di= -a ; di <= a; di ++) 



for(dj= - a; dj <= a ; dj ++) 

template [i+di,j+dj] = OCCUPIED 
end_f unction 

We have used the following parameters : NX = NY = 32, 
NPOINTS = 6, size = 3, and calculated the mean value 
of the Minkowski functionals over 1000 distributions. 
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